Module 1 Python
In [1]:
# Import Packages
from SimPEG import *
from IPython.html.widgets import interactive
import numpy as np
%pylab inline
In [2]:
# Begin by creating a ficticious set of model data
N=1000 # Set number of model parameters
x=np.linspace(0,1,N+1) # Define 1D domain on nodes
xc = 0.5*(x[1:] + x[0:-1]) # Define 1D domain on cell centers
m = np.zeros(N) # preallocate m array
# Define Gaussian function:
gauss = lambda x, a, mean, SD: a*np.exp(-(x-mean)**2./(2.*SD**2.)) # create a Gaussian function
# Choose parameters for Gaussian, evaluate, and store in an array, f.
SD=6
mean=0
a=1/(SD*sqrt(2*pi))
x2=np.linspace(-20,20,0.2*N)
f = gauss(x2,a, mean,SD)
# Define a boxcar function:
box = 0.4*np.ones(0.2*N)
# Insert the Gaussian and Boxcar into m:
m[0.2*N:0.4*N]=box
m[0.6*N:0.8*N]=10*f
# Plot
pylab.plot(xc,m)
#pylab.plot(f)
print size(f)
pylab.xlabel('x')
pylab.ylabel('m(x)')
pylab.title('Model, $m(x)$')
print shape(xc)
In [3]:
# Make the set of kernel functions
g = lambda x, i, p, q: np.exp(-p*i*x)*np.cos(2*np.pi*q*i*x) # create an anonymous function as immediately above
p = 0.01 # Set values for p, q
q = 0.15
M = 20 # specify number of output data
Gn = np.zeros((N+1,M))
for i in range(M):
f = g(x,i,p,q)
Gn[:,i] = f
# Plot
pylab.plot(x,Gn)
pylab.xlabel('x')
pylab.ylabel('g(x)')
pylab.title('Kernel functions, $g(x)$')
Out[3]:
In [4]:
# Make Averaging Matrix
n=N+1 # Define n as the N+1 dimension of the matrix
w=n-1 # Define w and the n-1 dimentsion of the matix
s = (N,n) # Store matrix values
Av = np.zeros(s) # Create a matrix of zeros of the correct dimensions
# and fill in with elements usin the loop below (note the 1/2 is included in here).
for i in range(N):
j=i
k=i+1
Av[i,j] = 0.5
Av[i,k] = 0.5
print shape(Av)
# Plot
pylab.title('Subset of averaging matrix')
pyplot.imshow(Av[:5,:6], cmap='Greys', interpolation='nearest')
print(Av[:5,:6])
In [5]:
# make the Volume, "delta x" array
Deltax = diff(x)
V = np.diag(Deltax) # create diagonal matrix
print shape(V)
print np.sum(Deltax)
print shape(x)
In [6]:
G=np.dot(np.transpose(np.dot(Av, Gn)), V)
d = np.dot(np.dot(np.transpose(np.dot(Av, Gn)), V),m) # Create synthetic data
xd=np.linspace(0,1,len(d)) # make x array for data
# Plot
pylab.plot(xd,d)
pylab.xlabel('')
pylab.ylabel('d')
pylab.title('Synthetic Data $d$')
Out[6]:
In [7]:
# Add noise to our synthetic data
noise = lambda fl, length, data, s: fl + randn(length)*data*s # introduce noise using a floor (f), length (l) and scaling factor (s)
r = noise(0, len(d), d, 0.04)
dobs= d + r
#pylab.plot(xd,r)
pylab.plot(xd,d)
pylab.plot(xd,dobs)
pylab.xlabel('')
pylab.ylabel('d')
pylab.title('Synthetic Data and noise, $d$')
Out[7]:
In [8]:
# Calculate the data misfit, phi_d #
####################################
invvar = lambda floor, percent, data: 1./(floor + percent*np.abs(data))
L2 = lambda M, w: dot(dot(w.T,M.T),dot(M,w)) # Anonymous function for 2-Norm
# assign a floor
flr = 0.015*dot(d.T,d)**0.5/len(d)
# Make Wd
eps = invvar(flr,0.02,dobs) # define epsilon and Wd
Wd=np.diag(eps)
# Take the 2-norm
phi_d = L2(Wd, np.dot(G,m)-dobs)
print phi_d
print flr
In [9]:
# Reference Model#
##################
## A constant reference model ##
mref = np.mean(m)*np.ones(N)
# Plot
pylab.plot(xc,mref)
pylab.xlabel('x')
pylab.ylabel('m(x)')
pylab.title('Reference Model, $m(x)$')
Out[9]:
In [10]:
# Domains #
############################### # Set up domains:
l = np.power(x[1:len(x)]-x[0:len(x)-1],0.5) # vector of distances between nodes, l
midx=np.dot(Av,x) # vector of midpoints, midx
h = np.power(midx[1:len(midx)]-x[0:len(midx)-1], 0.5) # Calculate distances between midpoints & take sqr roots of each value in vector h
In [11]:
# Calculate the model regularization, phi_m #
##############################################
# Create Ws, Wx matrices
Ws=np.diag(l) # put length vector into a diagonal matrix, Ws
Wx = np.zeros((len(m), len(m))) # preaollocate array and enter values into matrix, Wx
for i in range(len(h)):
j=i
k=i+1
Wx[i,j] = h[i]
Wx[i,k] = h[i]
alpha_s=0.1 # Set alpha_s and alpha_x values
alpha_x=0.15
# build Wm #
Wm=np.concatenate((alpha_s*Ws, alpha_x*Wx), axis=0)
beta = 1000000
# Set beta value
# Get recovered model
mrecovered = lambda G,Wd,Wm,data,mref, beta: dot(np.linalg.inv(dot(dot(G.T,Wd.T),dot(Wd,G)) + beta*dot(Wm.T,Wm)),dot(dot(G.T,Wd),dot(Wd,dobs)) + beta*dot(dot(Wm.T,Wm),mref))
# mrec = dot(np.linalg.inv(dot(dot(G.T,Wd.T),dot(Wd,G)) + beta*dot(Wm.T,Wm)),dot(dot(G.T,Wd),dot(Wd,dobs)) + beta*dot(dot(Wm.T,Wm),mref))
mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
# Get residual of mref and mrec
phi_d2 = L2(Wd,np.dot(G,mrec)-dobs)
phi_m2 = L2(Wm, mref-mrec)
print phi_d2
print phi_m2
# Plot
pylab.plot(xc,m)
pylab.plot(xc,mrec)
pylab.xlabel('x')
pylab.ylabel('m(x)')
pylab.title('Model, $m(x)$')
Out[11]:
In [58]:
# Get predicted data
dpred = np.dot(G,mrec)
#pylab.plot(xd,r)
pylab.plot(xd,dpred,'o')
pylab.plot(xd,dobs)
pylab.plot(xd,d,'x')
pylab.xlabel('')
pylab.ylabel('d')
pylab.title('Predicted and Observed Data')
Out[58]:
The $\beta$ Cooling method
mar 25
In [12]:
print phi_d
In [ ]:
# Make the Tikhonov Curve
beta = 100000 # Set beta value
itns = 10 # set maximum number of iterations
phdTik=np.zeros(itns) # preallocate arrays
phmTik=np.zeros(itns)
foo = np.zeros(itns)
b=np.zeros(itns)
for i in range(itns):
mrec = dot(np.linalg.inv(dot(dot(G.T,Wd.T),dot(Wd,G)) + beta*dot(Wm.T,Wm)),dot(dot(G.T,Wd),dot(Wd,dobs)) + beta*dot(dot(Wm.T,Wm),mref))
b[i] = beta
phdTik[i] = L2(Wd, np.dot(G,mrec)-dobs)
phmTik[i] = L2(Wm, mref-mrec)
foo[i] = np.absolute(phdTik[i]-phdstar)
beta=beta*0.5
pylab.plot(phmTik,phdTik,'o')
pylab.xlabel('$\phi_m$')
pylab.ylabel('$\phi_d$')
pylab.title('Tikhonov Curve')
This cell defines a function to find beta
In [76]:
def bisect_beta(tol, phid, phi_dstar, G,Wd,Wm,dobs,mref,beta):
i=1
while np.absolute(phid - phi_dstar)>tol:
mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
phid = L2(Wd,np.dot(G,mrec)-dobs)
print("Iteration = ",j, "beta=", beta, "phid at start=",phid)
if phid>phi_dstar:
beta = beta_0*0.5**i
print("beta=", beta)
if phid<phi_dstar:
beta_abv = beta_0*0.5**(i-2)
beta_bel = beta_0*0.5**(i-1)
beta_mid = 0.5*(beta_abv+beta_bel)
mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
phidmid = L2(Wd,np.dot(G,mmid)-dobs)
#for k in range(itns):
while np.absolute(phidmid - phi_dstar)>tol:
if phidmid<phi_dstar:
beta_abv = beta_0*0.5**(i-2)
beta_bel = beta_mid
beta_mid = 0.5*(beta_abv+beta_bel)
beta=beta_mid
mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
phidmid = L2(Wd,np.dot(G,mmid)-dobs)
if phidmid>phi_dstar:
beta_abv = beta_mid
beta_bel = beta_bel
beta_mid = 0.5*(beta_abv+beta_bel)
beta=beta_mid
mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
phidmid = L2(Wd,np.dot(G,mmid)-dobs)
if phidmid<phi_dstar:
beta_abv = beta_abv
beta_bel = beta_mid
beta_mid = 0.5*(beta_abv+beta_bel)
beta=beta_mid
mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
phidmid = L2(Wd,np.dot(G,mmid)-dobs)
phid=phidmid
mrec=mmid
return beta, phid, mrec
i=i+1
In [78]:
tol=0.01
phi
bisect_beta(tol, phid, phi_dstar, G,Wd,Wm,dobs,mref,beta)
In [ ]:
**This cell does the same as the function**
In [72]:
beta_0=1000000
beta=beta_0
tol=0.00001
i=1
itns=10
mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
phid = L2(Wd,np.dot(G,mrec)-dobs)
phi_dstar=20
while np.absolute(phid - phi_dstar)>tol:
j=i
mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
phid = L2(Wd,np.dot(G,mrec)-dobs)
print("Iteration = ",j, "beta=", beta, "phid at start=",phid)
if phid>phi_dstar:
beta = beta_0*0.5**i
print("beta=", beta)
if phid<phi_dstar:
beta_abv = beta_0*0.5**(i-2)
beta_bel = beta_0*0.5**(i-1)
beta_mid = 0.5*(beta_abv+beta_bel)
mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
phidmid = L2(Wd,np.dot(G,mmid)-dobs)
#for k in range(itns):
while np.absolute(phidmid - phi_dstar)>tol:
if phidmid<phi_dstar:
beta_abv = beta_0*0.5**(i-2)
beta_bel = beta_mid
beta_mid = 0.5*(beta_abv+beta_bel)
beta=beta_mid
mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
phidmid = L2(Wd,np.dot(G,mmid)-dobs)
if phidmid>phi_dstar:
beta_abv = beta_mid
beta_bel = beta_bel
beta_mid = 0.5*(beta_abv+beta_bel)
beta=beta_mid
mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
phidmid = L2(Wd,np.dot(G,mmid)-dobs)
if phidmid<phi_dstar:
beta_abv = beta_abv
beta_bel = beta_mid
beta_mid = 0.5*(beta_abv+beta_bel)
beta=beta_mid
mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
phidmid = L2(Wd,np.dot(G,mmid)-dobs)
phid=phidmid
mrec=mmid
i=i+1
print("At end beta=", beta,"phid=", phid)
This is the recovered model for the given beta
In [66]:
beta= 32462.310791015625
mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
phid = L2(Wd,np.dot(G,mrec)-dobs)
print phid
# Plot
pylab.plot(xc,m)
pylab.plot(xc,mrec)
pylab.xlabel('x')
pylab.ylabel('m(x)')
pylab.title('Model, $m(x)$')
Out[66]:
This is a version that shows what is going on inside ivery iteration
In [65]:
beta_0=100000
beta=beta_0
tol=0.1
i=1
itns=10
mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
phid = L2(Wd,np.dot(G,mrec)-dobs)
phi_dstar=20
while np.absolute(phid - phi_dstar)>tol:
j=i
mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
phid = L2(Wd,np.dot(G,mrec)-dobs)
print("Iteration = ",j, "beta=", beta, "phid at start=",phid)
if phid>phi_dstar:
beta = beta_0*0.5**i
print("beta=", beta)
if phid<phi_dstar:
beta_abv = beta_0*0.5**(i-2)
beta_bel = beta_0*0.5**(i-1)
beta_mid = 0.5*(beta_abv+beta_bel)
mabv = mrecovered(G,Wd,Wm,dobs,mref,beta_abv)
mbel = mrecovered(G,Wd,Wm,dobs,mref,beta_bel)
mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
phidabv = L2(Wd,np.dot(G,mabv)-dobs)
phidbel = L2(Wd,np.dot(G,mbel)-dobs)
phidmid = L2(Wd,np.dot(G,mmid)-dobs)
print("Itn # is:",j, "below threshold", "beta above=", beta_abv, "beta below=", beta_bel, "beta mid=",beta_mid)
print("Itn # is:",j, "below threshold", "phi above=", phidabv, "phi below=", phidbel, "phi mid=", phidmid)
for k in range(itns):
#while np.absolute(phidmid - phi_dstar)>tol:
if phidmid<phi_dstar:
beta_abv = beta_0*0.5**(i-2)
beta_bel = beta_mid
beta_mid = 0.5*(beta_abv+beta_bel)
beta=beta_mid
mabv = mrecovered(G,Wd,Wm,dobs,mref,beta_abv)
mbel = mrecovered(G,Wd,Wm,dobs,mref,beta_bel)
mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
phidabv = L2(Wd,np.dot(G,mabv)-dobs)
phidbel = L2(Wd,np.dot(G,mbel)-dobs)
phidmid = L2(Wd,np.dot(G,mmid)-dobs)
print("Itn # is:",j, "in R1", "beta above=", beta_abv, "beta below=", beta_bel, "beta mid=",beta_mid)
print("Itn # is:",j, "in R1", "phi above=", phidabv, "phi below=", phidbel, "phi mid=", phidmid)
if phidmid>phi_dstar:
beta_abv = beta_mid
beta_bel = beta_bel
beta_mid = 0.5*(beta_abv+beta_bel)
beta=beta_mid
mabv = mrecovered(G,Wd,Wm,dobs,mref,beta_abv)
mbel = mrecovered(G,Wd,Wm,dobs,mref,beta_bel)
mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
phidabv = L2(Wd,np.dot(G,mabv)-dobs)
phidbel = L2(Wd,np.dot(G,mbel)-dobs)
phidmid = L2(Wd,np.dot(G,mmid)-dobs)
print("Itn # is:",j, "in R2", "beta above=", beta_abv, "beta below=", beta_bel, "beta mid=",beta_mid)
print("Itn # is:",j, "in R2", "phi above=", phidabv, "phi below=", phidbel, "phi mid=", phidmid)
if phidmid<phi_dstar:
beta_abv = beta_abv
beta_bel = beta_mid
beta_mid = 0.5*(beta_abv+beta_bel)
beta=beta_mid
mabv = mrecovered(G,Wd,Wm,dobs,mref,beta_abv)
mbel = mrecovered(G,Wd,Wm,dobs,mref,beta_bel)
mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
phidabv = L2(Wd,np.dot(G,mabv)-dobs)
phidbel = L2(Wd,np.dot(G,mbel)-dobs)
phidmid = L2(Wd,np.dot(G,mmid)-dobs)
print("Itn # is:",j, "in R3", "beta above=", beta_abv, "beta below=", beta_bel, "beta mid=",beta_mid)
print("Itn # is:",j, "in R3", "phi above=", phidabv, "phi below=", phidbel, "phi mid=", phidmid)
phid=phidmid
mrec=mmid
i=i+1
print("At end beta=", beta,"phid=", phid)